##Robert Schub and Allen Schmaltz
##Proactive Reputation Building
##This file contains matching tests of Clare and Danilovic (2010) and work on general entry deterrence


library(foreign)
library(xtable)
library(car)
library(robustbase)
library(mgcv)
library(MASS)
library(Zelig)
library(survey)
library(lme4)
library(VGAM)
library(mlogit)
library(sandwich)
library(gee)
library(Design)
library(lpSolve)
library(WhatIf)
library(lmtest)
library(zoo)
library(MatchIt)
library(tcltk)
library(cem)



cla <- read.dta(file="clare.dta")
head(cla)
dim(cla)

##########################################################
##### RUN ORIGINAL MODEL 1.2 W/ NO PRUNING ###############
##########################################################

cl.orig<-cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])

model1.2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = cla)

#Get clustered standard errors
mat2.2<-estfun(model1.2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.orig$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.orig$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.orig<-df2.2*sandwich(model1.2,meat=crossprod(u2.2)/N2.2)
(vcov.orig)
ses.orig<-coeftest(model1.2,vcov.orig)
ses.orig

##RESULTS##
table.orig<-ses.orig[,1:2]
xtable(table.orig,digits=3)


##########################################################
##### 2x2 MATCHING SCHEME: 3 BOXES AS CONTROL ############
##### 1 AS TREATMENT:::USE CEM MATCHING       ############
##########################################################

## Grab revelant data for Model 1.2
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)

#Create Treatment Variable 
cl1.2$treat<-cl1.2$weakmr_decline9
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1>0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1>0]<-1
hist(cl1.2$treat)

head(cl1.2)
cem.2<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad"))
cem.2
## Keep 5105 Control and 1222 Treated with L1 of 0.600 ##
(5105+1222)/12192


#Add variable noting whether that observation was matched and what weight to put on it
head(cl1.2)
cl1.2$IN <- cem.2$matched
weights<-cem.2$w
cl1.2<-cbind(cl1.2,weights)


#Create dataset with only matched observation
cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]
dim(cl1.2IN)

#Run parametric model on the matched dataset with proper weights
cem.W<-zelig(formula=cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights)
summary(cem.W)

#Get clustered standard errors
mat2.2<-estfun(cem.W)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL2.2<-df2.2*sandwich(cem.W,meat=crossprod(u2.2)/N2.2)
(vcovCL2.2)
ses2.2<-coeftest(cem.W,vcovCL2.2)
ses2.2


##RESULTS##
table2<-ses2.2[,1:2]
xtable(table2,digits=3)
#Interaction terms has almost no effect and is insignificant#


##########################################################
##### DEMONSTRATION OF BALANCE: PRE & POST MATCHING ######
##########################################################

#Use matchit to generate balance statistics
auto.match<-matchit(formula=treat ~ weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl1.2, method="cem")
summary(auto.match)

#Pulling key balance measures
pre.cem<-summary(auto.match)$sum.all
post.cem<-summary(auto.match)$sum.matched

pre.stdc<-(pre.cem[,1]-pre.cem[,2])/pre.cem[,3]
post.stdc<-(post.cem[,1]-post.cem[,2])/post.cem[,3]
pre.stdc

rows<-c("Distance","TargetRep","RelCap","JointDem","PolRel","PceYrs","PceYrs2","PceYrs3")

#Graph of balance improvement
pdf(file="balance.pdf")
plot(x = pre.stdc, y = 1:nrow(pre.cem),
     pch = 19, xlab = "Standardized Imbalance: CEM", axes = FALSE, 
     xlim = c(-0.2,0.2), ylab = "")
points(x = post.stdc, y = 1:nrow(post.cem),
     pch = 17)
abline(v = 0, lty = "dashed")
axis(1)
axis(2, at = 1:8, labels=rows,las = 1, cex.axis = .8)
legend("topright", legend = c("Pre-match:  L1=0.779","Post-match: L1=0.600"), pch = c(19,17), cex = .8)
dev.off()




##########################################################
##### FIRST DIFFERENCES: ORIGINAL VS. MATCHED PLOT #######
##########################################################

## For this analysis I will first grab the medians for all control variables from the MATCHED data, then run first differences moving rivals, weak rep, and the interaction term from all 0s to the scenario of rivals=XX, weak rep=XX, and interaction term is the product of these two

## Then, I use these same medians and estimate the first difference using the original model used by Clare/Danilovic. The two first differences (density plots) are then plotted against one another


#Creating covariates 
cov<-apply(X=cl1.2IN,MARGIN=2,FUN=median)
cov<-c(1,cov[2:11])

#Creating control covariates
cov.con<-cov
cov.con[2]<-0
cov.con[3]<-0
cov.con[4]<-cov.con[2]*cov.con[3]

#Creating treated covariates [Maxing out weak and rivals]
cov.treat<-cov
cov.treat[2]<-6
cov.treat[3]<-0.9
cov.treat[4]<-cov.treat[2]*cov.treat[3]

#MATCHED#
#Grab coefficients from MATHCED model and draw beta from the MVN
beta<-table2[,1]
beta.draw<-mvrnorm(100000,mu=beta,Sigma=vcovCL2.2)

#Create expected values for treat and control and subtract
est.con<-pnorm(cov.con%*%t(beta.draw))
est.treat<-pnorm(cov.treat%*%t(beta.draw))
first.diff.match<-(est.treat-est.con)
mean(first.diff.match)

plot(density(est.con))
lines(density(est.treat),lty=2)

#ORIGINAL#
#Grab coefficients from ORIGINAL model and draw beta from the MVN
beta.orig<-table.orig[,1]
beta.orig.draw<-mvrnorm(100000,mu=beta.orig,Sigma=vcov.orig)

#Create expected values for treat and control and subtract
est.con.orig<-pnorm(cov.con%*%t(beta.orig.draw))
est.treat.orig<-pnorm(cov.treat%*%t(beta.orig.draw))
first.diff.orig<-(est.treat.orig-est.con.orig)
mean(first.diff.orig)

#GRAPH First Diffs
pdf(file="firstdiffs.pdf")
plot(density(first.diff.orig),xlim=c(-0.05,0.1),ylim=c(0,25), main="Increased Probability of Dispute Initiation*",xlab="*Moving from Control to Maximum Treament Dosage")
lines(density(first.diff.match),lty=2)
legend("topright", legend = c("Naive Data","Matched Data"), lty=c(1,2))
abline(v=0,lty=3)
dev.off()


##########################################################
##### FIGURE 1 OF ORIGINAL VS. MATCHED PLOTS #############
##########################################################

################ FIGURE 1 MATCHED ######################
#cem.W
#vcovCL2.2
cem.W
library(MASS)
#beta.draws <- mvrnorm(10000, mu = as.matrix(table1Model2$coefficients)[,1], Sigma = summary(table1Model2)$cov.unscaled)
beta.draws <- mvrnorm(10000, mu = as.matrix(cem.W$coefficients)[,1], Sigma = vcovCL2.2)

#NEED COVARATIES TO BE MEDIANS FROM MATCHED DATA: NOT SAME AS ORIGINAL AUTHORS#
covs<-apply(X=cl1.2IN,MARGIN=2,FUN=median)
covs

highrisk3 <- rep(0, 11)
#highrisk3[2:11] <- highrisk
highrisk3[1] <- 1
#values the original authors used in Stata:
highrisk3[5] <- covs[5]
highrisk3[6] <- covs[6]
highrisk3[7] <- covs[7]
highrisk3[8] <- 1  #scalar h_PolRel=1
highrisk3[9] <- covs[9]
highrisk3[10] <- covs[10]
highrisk3[11] <- covs[11]

highrisk3 <- as.matrix(highrisk3)
dim(t(beta.draws))
#test
#p.ests2 <- pnorm(t(highrisk3)%*%t(beta.draws))
#

NUMBERrivals <- seq(0, 6, 2)
yrsWeak <- 1:10
pEsts0 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts2 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts4 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts6 <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)

for(NUMBERrivalsIndex in 1:length(NUMBERrivals)){
  for(j in 1:length(yrsWeak)){
    highriskweak <- highrisk3
    highriskweak[2] <- NUMBERrivals[NUMBERrivalsIndex]
    highriskweak[3] <- 0.9^j
    highriskweak[4] <- NUMBERrivals[NUMBERrivalsIndex] * 0.9^j
    #highriskweak[3] <- yrsWeak[j]
    #highriskweak[4] <- NUMBERrivals * yrsWeak[j]
    if (NUMBERrivals[NUMBERrivalsIndex]==0) {
    	  pEsts0[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	} 
    if (NUMBERrivals[NUMBERrivalsIndex]==2) {
    	  pEsts2[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==4) {
    	  pEsts4[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==6) {
    	  pEsts6[,j] <- pnorm(t(highriskweak)%*%t(beta.draws))
    	}       	
  }
}

#sanity check
#class(highriskweak)
#highriskweak<-as.numeric(highriskweak)
#highriskweak[2]<-2
#highriskweak[3]<-0.4
#highriskweak[4]<-highriskweak[2]*highriskweak[3]
#cof<-colMeans(beta.draws)
#mean(pnorm(highriskweak%*%t(beta.draws)))
#colMeans(pEsts2)
#

pdf(file="fig1mat.pdf")
plot(yrsWeak, apply(pEsts0,2,mean), ylim = c(0,0.2), lty=1, ylab="Probability of Initiation", xlab="Years Since Challenger Backed Down", type="l",main="Predicted Probability of Initiation: Matched Data")
segments(x0 = yrsWeak+0.1, x1 = yrsWeak+0.1,
y0 = apply(pEsts0, 2, quantile, .025),
y1 = apply(pEsts0, 2, quantile, .975), lty=1)
lines(x = yrsWeak, y = apply(pEsts2,2,mean),lty=2)
segments(x0 = yrsWeak+0.05, x1 = yrsWeak+0.05,
y0 = apply(pEsts2, 2, quantile, .025),
y1 = apply(pEsts2, 2, quantile, .975), lty=2)
lines(x = yrsWeak, y = apply(pEsts4,2,mean), lty=3)
segments(x0 = yrsWeak+0.15, x1 = yrsWeak+0.15,
y0 = apply(pEsts4, 2, quantile, .025),
y1 = apply(pEsts4, 2, quantile, .975), lty=3)
lines(x = yrsWeak, y = apply(pEsts6,2,mean),lty=4)
segments(x0 = yrsWeak, x1 = yrsWeak,
y0 = apply(pEsts6, 2, quantile, .025),
y1 = apply(pEsts6, 2, quantile, .975),lty=4)
legend("topright", c("0 Other Rivals","2 Other Rivals", "4 Other Rivals", "6 Other Rivals"), lwd=c(1, 1, 1, 1), lty=c(1, 2, 3, 4))
dev.off()

################ FIGURE 1 ORIGINAL ######################

cl.orig<-cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])

model1.2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = cla)

#Get clustered standard errors
mat2.2<-estfun(model1.2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.orig$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.orig$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.orig<-df2.2*sandwich(model1.2,meat=crossprod(u2.2)/N2.2)
(vcov.orig)
ses.orig<-coeftest(model1.2,vcov.orig)
ses.orig

library(MASS)
#beta.draws <- mvrnorm(10000, mu = as.matrix(table1Model2$coefficients)[,1], Sigma = summary(table1Model2)$cov.unscaled)
beta.draws.o <- mvrnorm(10000, mu = as.matrix(model1.2$coefficients)[,1], Sigma = vcov.orig)

#NEED COVARATIES TO BE MEDIANS FROM MATCHED DATA: NOT SAME AS ORIGINAL AUTHORS#
covs<-apply(X=cl1.2IN,MARGIN=2,FUN=median)
covs

highrisk3 <- rep(0, 11)
#highrisk3[2:11] <- highrisk
highrisk3[1] <- 1
#values the original authors used in Stata:
highrisk3[5] <- covs[5]
highrisk3[6] <- covs[6]
highrisk3[7] <- covs[7]
highrisk3[8] <- 1  #scalar h_PolRel=1
highrisk3[9] <- covs[9]
highrisk3[10] <- covs[10]
highrisk3[11] <- covs[11]

highrisk3 <- as.matrix(highrisk3)
dim(t(beta.draws.o))
#test
#p.ests2 <- pnorm(t(highrisk3)%*%t(beta.draws))
#

NUMBERrivals.o <- seq(0, 6, 2)
yrsWeak <- 1:10
pEsts0.o <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts2.o <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts4.o <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)
pEsts6.o <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)

for(NUMBERrivalsIndex in 1:length(NUMBERrivals.o)){
  for(j in 1:length(yrsWeak)){
    highriskweak <- highrisk3
    highriskweak[2] <- NUMBERrivals[NUMBERrivalsIndex]
    highriskweak[3] <- 0.9^j
    highriskweak[4] <- NUMBERrivals[NUMBERrivalsIndex] * 0.9^j
    #highriskweak[3] <- yrsWeak[j]
    #highriskweak[4] <- NUMBERrivals * yrsWeak[j]
    if (NUMBERrivals[NUMBERrivalsIndex]==0) {
    	  pEsts0.o[,j] <- pnorm(t(highriskweak)%*%t(beta.draws.o))
    	} 
    if (NUMBERrivals[NUMBERrivalsIndex]==2) {
    	  pEsts2.o[,j] <- pnorm(t(highriskweak)%*%t(beta.draws.o))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==4) {
    	  pEsts4.o[,j] <- pnorm(t(highriskweak)%*%t(beta.draws.o))
    	}
    if (NUMBERrivals[NUMBERrivalsIndex]==6) {
    	  pEsts6.o[,j] <- pnorm(t(highriskweak)%*%t(beta.draws.o))
    	}       	
  }
}


pdf(file="fig1orig.pdf")
plot(yrsWeak, apply(pEsts0.o,2,mean), ylim = c(0,0.2), type="l", ylab="Probability of Initiation", xlab="Years Since Challenger Backed Down",lty=1,main="Predicted Probability of Initiation: Naive Data")
segments(x0 = yrsWeak+0.1, x1 = yrsWeak+0.1,
y0 = apply(pEsts0.o, 2, quantile, .025),
y1 = apply(pEsts0.o, 2, quantile, .975), lty=1)
lines(x = yrsWeak, y = apply(pEsts2.o,2,mean), lty=2)
segments(x0 = yrsWeak+0.05, x1 = yrsWeak+0.05,
y0 = apply(pEsts2.o, 2, quantile, .025),
y1 = apply(pEsts2.o, 2, quantile, .975), lty=2)
lines(x = yrsWeak, y = apply(pEsts4.o,2,mean), lty=3)
segments(x0 = yrsWeak+0.15, x1 = yrsWeak+0.15,
y0 = apply(pEsts4.o, 2, quantile, .025),
y1 = apply(pEsts4.o, 2, quantile, .975), lty=3)
lines(x = yrsWeak, y = apply(pEsts6.o,2,mean),lty=4)
segments(x0 = yrsWeak, x1 = yrsWeak,
y0 = apply(pEsts6.o, 2, quantile, .025),
y1 = apply(pEsts6.o, 2, quantile, .975),lty=4)
legend("topright", c("0 Other Rivals","2 Other Rivals", "4 Other Rivals", "6 Other Rivals"),lwd=c(1, 1, 1, 1), lty=c(1, 2, 3, 4))
dev.off()




##########################################################
##### ALTERNATIVE PRUNING/MODEL DEPENDENCE CHECKS ########
##########################################################

##########################################################
##### (1) 3x3 MATCHING TO MATCH ACROSS DOSAGE LEVELS #####
##########################################################
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)

#Create treatment variable
cl1.2$weak<-cl1.2$weakmr_decline9
cl1.2$weak[cl1.2$weakmr_decline9==0]<-0
cl1.2$weak[cl1.2$weakmr_decline9 > 0 & cl1.2$weakmr_decline9 < 0.7]<-1
cl1.2$weak[cl1.2$weakmr_decline9 > .7]<-2
hist(cl1.2$weak)

cl1.2$rival<-cl1.2$rivalry_count1
cl1.2$rival[cl1.2$rivalry_count1==0]<-0
cl1.2$rival[cl1.2$rivalry_count1 > 0 & cl1.2$rivalry_count1 < 3]<-1
cl1.2$rival[cl1.2$rivalry_count1 >2]<-2
hist(cl1.2$rival)

cl1.2$treat<-cl1.2$weak
cl1.2$treat[cl1.2$weak==0 | cl1.2$rival==0]<-0
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==1]<-1
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==1]<-2
cl1.2$treat[cl1.2$weak==1 & cl1.2$rival==2]<-3
cl1.2$treat[cl1.2$weak==2 & cl1.2$rival==2]<-4

#Match across control and all four treatments
cem.mat3<-cem(treatment="treat", data=cl1.2, drop=c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","dyad","weak","rival"))
cem.mat3

##LEFT WITH 502C, 54T1, 33T2, 43T3, 44T4##
##L1 = 0.636

#Creating subsetted data w/ only matched observations
cl1.2$IN <- cem.mat3$matched
weights3<-cem.mat3$w
cl1.2<-cbind(cl1.2,weights3)
cl1.2IN <- cl1.2[cl1.2$IN==TRUE,]

#Run parametric model with only matched data and proper weights
cem.3<-zelig(formula=cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + relcap + jdem + pceyrs + pceyrs2 + pceyrs3, data=cl1.2IN, model="probit", weights=cl1.2IN$weights3)
summary(cem.3)

##sandwich chokes bc no variation in pol_rel and weakmr2 so drop both from regression
mat2.2<-estfun(cem.3)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl1.2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl1.2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcovCL3<-df2.2*sandwich(cem.3,meat=crossprod(u2.2)/N2.2)
(vcovCL3)
ses3<-coeftest(cem.3,vcovCL3)
ses3

table3<-ses3[,1:2]
xtable(table3, caption="3x3 CEM",digits=3)


##########################################################
##### (2) CONVEX HULL PRUNING  ###########################
##########################################################
## Grab revelant data for Model 1.2
head(cla)
cl1.2<-na.omit(cla[c("cwinit","rivalry_count1","weakmr_decline9","nriv_weakdec9mr","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])
dim(cl1.2)

#Resetting rivalry count 1 and the interaction term
cl1.2$rivalry_count1_orig<-cl1.2$rivalry_count1
cl1.2$rivalry_count1<-cl1.2$rivalry_count1_orig-1

cl1.2$nriv_weakdec9mr_orig<-cl1.2$nriv_weakdec9mr
cl1.2$nriv_weakdec9mr<-cl1.2$rivalry_count1 * cl1.2$weakmr_decline9

#Drop original terms here
dim(cl1.2)
cl1.2<-cl1.2[,-14]
cl1.2<-cl1.2[,-13]
head(cl1.2)

#Create Treatment Variable 
cl1.2$treat<-cl1.2$weakmr_decline9
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9==0 & cl1.2$rivalry_count1>0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1==0]<-0
cl1.2$treat[cl1.2$weakmr_decline9>0 & cl1.2$rivalry_count1>0]<-1

## Convex Hull: Flip "treat" from 0 to 1 and 1 to 0 ##




## Hull 1: Flip 0 to 1 ##

##NEED TO RUN THIS LINE OF CODE!!!. MARKED OFF BECAUSE IT REQUIRES EXCESS PROCESSING TIME##

#a<-whatif(~weakmr2_decline9+relcap+jdem+pol_rel+pceyrs, data=cl1.2[cl1.2$treat==1,], cfact=cl1.2[cl1.2$treat==0,]) 
names(a)
a$in.hull
a.in<-(a[a$in.hull==TRUE])
a.out<-a[a$in.hull==FALSE]
length(a.in)
length(a.out)
## 8841 out of 10655 have a counterfactual in the hull ##

## Hull 2: Flip 1 to 0 ##

##NEED TO RUN THIS LINE OF CODE!!!. MARKED OFF BECAUSE IT REQUIRES EXCESS PROCESSING TIME##

#b<-whatif(~weakmr2_decline9+relcap+jdem+pol_rel+pceyrs, data=cl1.2[cl1.2$treat==0,], cfact=cl1.2[cl1.2$treat==1,]) 
summary(b)
## 1488 out of 1537 have a counterfactual in the hull ##
b.in<-(b[b$in.hull==TRUE])
length(b.in)

#Dividing original relevant dataset by treatment condition#
data1.2.0<-cl1.2[cl1.2$treat==0,]
data1.2.1<-cl1.2[cl1.2$treat==1,]

#Removing observations not in convex hull 1#
bh<-b$in.hull
bh[bh==TRUE]<-1
turn.off1<-cbind(bh,data1.2.1)
turn.off1<-turn.off1[turn.off1[,1]==1,]
dim(turn.off1)

#Removing observations not in convex hull 2#
ah<-a$in.hull
ah[ah==TRUE]<-1
turn.off0<-cbind(ah,data1.2.0)
turn.off0<-turn.off0[turn.off0[,1]==1,]
dim(turn.off0)

#Reconstructing the complete dataset with only convex hull compliant observations#
data.1.2.pruned<-rbind(turn.off1[,-1],turn.off0[,-1])
dim(data.1.2.pruned)

#Rerun model with only the convex hull compliant dataset#
ch1.2 <- zelig(cwinit ~ rivalry_count1 + weakmr_decline9 + nriv_weakdec9mr + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, model = "probit", data = data.1.2.pruned)
summary(ch1.2)

#Get clustered standard errors
mat2.2<-estfun(ch1.2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,data.1.2.pruned$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(data.1.2.pruned$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.ch<-df2.2*sandwich(ch1.2,meat=crossprod(u2.2)/N2.2)
(vcov.ch)
ses.ch<-coeftest(ch1.2,vcov.ch)
ses.ch

##RESULTS##
table.ch<-ses.ch[,1:4]
xtable(table.ch,digits=3)

##EVEN WITH SUCH LIGHT PRUNING WE LOSE THE EFFECT--REDUCES SIZE w/ VARIANCE BASICALLY STABLE



################*************##############################
################*************##############################




##########################################################
#####            ENTRY DETERRENCE WORK        ############
##########################################################

##Pull all data for model 1.1 and 1.2 but now weak2 is the core treatment and will also run an interaction between weak2 and rival count2##

############
## TESTING WEAK REPUTATION OF TARGET:: NO INTERACTION #########
############
cl.wk1<-na.omit(cla[c("cwinit","rivalry_count2","weakmr_decline9","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])

#Reset rivalry
#Resetting rivalry count 2 
cl.wk1$rivalry_count2_orig<-cl.wk1$rivalry_count2
cl.wk1$rivalry_count2<-cl.wk1$rivalry_count2_orig-1
#Delete excess rival count column
cl.wk1<-cl.wk1[,-12]

#Run probit without matching
wk1<-zelig(cwinit~rivalry_count2+weakmr_decline9+weakmr2_decline9+relcap+jdem+pol_rel+pceyrs+pceyrs2+pceyrs3, data=cl.wk1, model="probit")
summary(wk1)

#Run sandwich
mat2.2<-estfun(wk1)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk1$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk1$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk1<-df2.2*sandwich(wk1,meat=crossprod(u2.2)/N2.2)
(vcov.wk1)
ses.wk1<-coeftest(wk1,vcov.wk1)
ses.wk1

table.wk1<-ses.wk1[,1:2]
xtable(table.wk1,digits=3,caption="naive weak rep (no interaction)")
## FIND WEAK REP IS SIGNIFICANT TO GETTING TARGETED (0.255 and SE 0.096)

##### NOW MATCH FOR 1.1 WEAK TARGET #####################

##NOW CEM FOR WEAK REP INTO TWO BUCKETS THEN RE-RUN#
head(cl.wk1)
#generate weak2 rep dummy#
cl.wk1$weak<-cl.wk1$weakmr2_decline9
cl.wk1$weak[cl.wk1$weakmr2_decline>0]<-1

cem.wk<-cem(treatment="weak", data=cl.wk1, drop=c("cwinit","weakmr2_decline9","dyad"))
cem.wk
## Keep 2117 Control and 1015 Treated with L1 of 0.676 ##

#Add variable noting whether that observation was matched and what weight to put on it
cl.wk1$IN <- cem.wk$matched
weight.wk<-cem.wk$w
cl.wk1<-cbind(cl.wk1,weight.wk)


#Create dataset with only matched observation
cl.wk1IN <- cl.wk1[cl.wk1$IN==TRUE,]
dim(cl.wk1IN)

#Run parametric model on the matched dataset with proper weights
cem.wk1.M<-zelig(formula=cwinit ~ rivalry_count2 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3, data=cl.wk1IN, model="probit", weights=cl.wk1IN$weight.wk)
summary(cem.wk1.M)

#Get clustered standard errors
mat2.2<-estfun(cem.wk1.M)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk1IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk1IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk1.M<-df2.2*sandwich(cem.wk1.M,meat=crossprod(u2.2)/N2.2)
(vcov.wk1.M)
ses.wk1.M<-coeftest(cem.wk1.M,vcov.wk1.M)
ses.wk1.M

table.wk1.M<-ses.wk1.M[,1:2]
xtable(table.wk1.M,digits=3,caption="weak2 rep matched has an effect")
##SOME EFFECT FOR WK2 POST-MATCHING (.217 and SE .130)


##########################################################
##### FIRST DIFFERENCES: MATCHED PLOT ####################
##########################################################

## For this analysis I will first grab the medians for all control variables from the MATCHED data, then run first differences moving weak2 rep from 0 to weak2 rep=0.9

#Creating covariates 
cov.wk<-apply(X=cl.wk1IN,MARGIN=2,FUN=median)
cov.wk<-c(1,cov.wk[2:10])
cov.wk
#Creating control covariates
cov.con.wk<-cov.wk
cov.con.wk[4]<-0

#Creating treated covariates [Maxing out weak and rivals]
cov.treat.wk<-cov.wk
cov.treat.wk[4]<-0.9

#Grab coefficients from MATHCED model and draw beta from the MVN
beta.wk<-table.wk1.M[,1]
beta.draw.wk<-mvrnorm(100000,mu=beta.wk,Sigma=vcov.wk1.M)

#Create expected values for treat and control and subtract
est.con.wk<-pnorm(cov.con.wk%*%t(beta.draw.wk))
est.treat.wk<-pnorm(cov.treat.wk%*%t(beta.draw.wk))
first.diff.match.wk<-(est.treat.wk-est.con.wk)
mean(first.diff.match.wk)

plot(density(est.con.wk))
lines(density(est.treat.wk),lty=2)

plot(density(first.diff.match.wk))

#GRAPH First Diffs
pdf(file="firstdiffsWEAK2.pdf")
plot(density(first.diff.match.wk),xlim=c(-0.05,0.1),ylim=c(0,25), main="Weak Target: Increased Probability of Dispute Initiation*",xlab="*Moving from Control to Maximum Treament Dosage")
abline(v=0,lty=3)
dev.off()


##########################################################
##### FIGURE 1 OF WEAK2 VARIED ACROSS LEVELS #############
##########################################################

################ FIGURE 1 MATCHED ######################
#cem.wk1.M
#vcov.wk1
cem.wk1.M
library(MASS)
#beta.draws <- mvrnorm(10000, mu = as.matrix(table1Model2$coefficients)[,1], Sigma = summary(table1Model2)$cov.unscaled)
beta.draws.wk <- mvrnorm(10000, mu = as.matrix(cem.wk1.M$coefficients)[,1], Sigma = vcov.wk1)

#NEED COVARATIES TO BE MEDIANS FROM MATCHED DATA: NOT SAME AS ORIGINAL AUTHORS#
cov.wk

covs.wk <- as.matrix(cov.wk)
dim(t(beta.draws.wk))
#test
#p.ests2 <- pnorm(t(highrisk3)%*%t(beta.draws))
#

NUMBERrivals <- seq(0, 6, 2)
yrsWeak <- 1:10
pEsts2.wk <- matrix(data = NA, ncol = length(yrsWeak), nrow=10000)

for(NUMBERrivalsIndex in 1:length(NUMBERrivals)){
  for(j in 1:length(yrsWeak)){
    covsweak <- covs.wk
    covsweak[2] <- NUMBERrivals[NUMBERrivalsIndex]
    covsweak[4] <- 0.9^j
    #highriskweak[3] <- yrsWeak[j]
    #highriskweak[4] <- NUMBERrivals * yrsWeak[j]
    if (NUMBERrivals[NUMBERrivalsIndex]==2) {
    	  pEsts2.wk[,j] <- pnorm(t(covsweak)%*%t(beta.draws.wk))
    	}       	
  }
}

#sanity check and creation of the 0 line (ie if not weak at all)
class(highriskweak)
covsweak<-as.numeric(covsweak)
covsweak[4]<-0
cof<-colMeans(beta.draws.wk)
control.wk<-mean(pnorm(covsweak%*%t(beta.draws.wk)))
colMeans(pEsts2.wk)

pdf(file="weaktarget.pdf")
plot(yrsWeak, apply(pEsts2.wk,2,mean), ylim = c(0.075,0.2), lty=1, ylab="Probability of Initiation", xlab="Years Since Target Backed Down", type="l",main="Predicted Probability of Initiation: Matched Data")
lines(x=yrsWeak,y=apply(pEsts2.wk,2,quantile,0.025), lty=2)
lines(x=yrsWeak,y=apply(pEsts2.wk,2,quantile,0.975), lty=2)
abline(h=control.wk,lty=3)
legend("topright", legend = c("Predicted Probability","95% Confidence Interval","Predicted Probability if Strong (Control Case)"), lty=c(1,2,4))
dev.off()


#############################################
## TESTING WEAK REPUTATION OF TARGET:: WITH INTERACTION ##############################################
##############
cl.wk2<-na.omit(cla[c("cwinit","rivalry_count2","weakmr_decline9","weakmr2_decline9","relcap","jdem","pol_rel","pceyrs","pceyrs2","pceyrs3","dyad")])

#Reset rivalry
#Resetting rivalry count 2 
cl.wk2$rivalry_count2_orig<-cl.wk2$rivalry_count2
cl.wk2$rivalry_count2<-cl.wk2$rivalry_count2_orig-1
#Delete excess rival count column
cl.wk2<-cl.wk2[,-12]

#Run naive probit
wk2<-zelig(cwinit~rivalry_count2+weakmr_decline9+weakmr2_decline9+relcap+jdem+pol_rel+pceyrs+pceyrs2+pceyrs3 + rivalry_count2*weakmr2_decline9, data=cl.wk1, model="probit")
summary(wk2)

#Run sandwich
mat2.2<-estfun(wk2)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk2$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk2$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk2<-df2.2*sandwich(wk2,meat=crossprod(u2.2)/N2.2)
(vcov.wk2)
ses.wk2<-coeftest(wk2,vcov.wk2)
ses.wk2

table.wk2<-ses.wk2[,1:2]
xtable(table.wk2,digits=3,caption="naive weak rep (with interaction)")
## FIND INTERACTION ISN'T SIGNIFICANT BUT WEAK REP LOWER ORDER IS SIGNIFICANT TO GETTING TARGETED (0.368 and SE 0.148)!

##### NOW MATCH FOR 1.2 WEAK TARGET #####################

##NOW CEM INTERACTION INTO TWO BUCKETS THEN RE-RUN#
head(cl.wk2)
#generate weak2 rep dummy#
cl.wk2$int<-cl.wk2$weakmr2_decline9
cl.wk2$int[cl.wk2$weakmr2_decline9==0 | cl.wk2$rivalry_count2==0]<-0
cl.wk2$int[cl.wk2$weakmr2_decline9>0 & cl.wk2$rivalry_count2>0]<-1

cem.wk2<-cem(treatment="int", data=cl.wk2, drop=c("cwinit","weakmr2_decline9","dyad","int","rivalry_count2"))
cem.wk2
## Keep 4813 Control and 1202 Treated with L1 of 0.593 ##

#Add variable noting whether that observation was matched and what weight to put on it
cl.wk2$IN <- cem.wk2$matched
weight.wk2<-cem.wk2$w
cl.wk2<-cbind(cl.wk2,weight.wk2)
dim(cl.wk2[cl.wk2$IN==TRUE,])

#Create dataset with only matched observation
cl.wk2IN <- cl.wk2[cl.wk2$IN==TRUE,]
dim(cl.wk2IN)

#Run parametric model on the matched dataset with proper weights
cem.wk2.M<-zelig(formula=cwinit ~ rivalry_count2 + weakmr_decline9 + weakmr2_decline9 + relcap + jdem + pol_rel + pceyrs + pceyrs2 + pceyrs3 + rivalry_count2*weakmr2_decline9, data=cl.wk2IN, model="probit", weights=cl.wk2IN$weight.wk2)
summary(cem.wk2.M)

#Get clustered standard errors
mat2.2<-estfun(cem.wk2.M)
mat2.2<-na.omit(mat2.2)
head(mat2.2)
dim(mat2.2)
N2.2<-nrow(mat2.2)
u2.2<-apply(mat2.2,2,function(x) tapply (x,cl.wk2IN$dyad,sum))
u2.2<-na.omit(u2.2)
clust2.2<-length(unique(cl.wk2IN$dyad))
a2.2<-clust2.2
b2.2<-nrow(mat2.2)
c2.2<-ncol(mat2.2)
df2.2<-((a2.2)/(a2.2-1))*((b2.2-1)/(b2.2-c2.2))
df2.2
vcov.wk2.M<-df2.2*sandwich(cem.wk2.M,meat=crossprod(u2.2)/N2.2)
(vcov.wk2.M)
ses.wk2.M<-coeftest(cem.wk2.M,vcov.wk2.M)
ses.wk2.M

table.wk2.M<-ses.wk2.M[,1:2]
xtable(table.wk2.M,digits=3,caption="weak2 rep matched has an effect BUT INT DOESN'T")
##NO EFFECT FOR INTERACTION POST-MATCHING BUT SOME FOR WEAK2 REP































